﻿using LxBasic;
using LxLESS;
using System.Diagnostics;
using System.Reflection.Emit;

namespace LxCore
{
    public abstract class SIMPLESolver
    {
        public class Set : IterateSolver.Set<Unit>
        {
            private int _NonOrthCorLoops = 0;
            public int NonOrthCorLoops
            {
                get => _NonOrthCorLoops;
                set => _NonOrthCorLoops = value >= 0 ? value : 0;
            }
            private int _PISOLoops = 0;
            public int PISOLoops
            {
                get => _PISOLoops;
                set => _PISOLoops = value >= 0 ? value : 0;
            }
            private int _AdditionalOuterLoops = 0;
            public int AdditionalOuterLoops
            {
                get => _AdditionalOuterLoops;
                set => _AdditionalOuterLoops = value >= 0 ? value : 0;
            }
            public override void CheckCreate()
            {
                base.CheckCreate();
                var list = new ConnectivitySolver().GetConnectivity<SIMPLEFields, EquationConnectBase>(f => f.SimpleGroup, c => null);//后面补充
                Solvers = new Unit[list.Count];
                for (int i1 = 0; i1 < Solvers.Length; i1++)
                {
                    var p = list[i1];
                    var s = new Unit();
                    s.CheckCreate(this, p.Equations.ToArray(), p.Connects.ToArray());
                    Solvers[i1] = s;
                }
            }
        }
        public class Unit : IterateSolver.Unit
        {
            public SIMPLEFields[] Equations;
            public EquationConnectBase[] Connects;
            public Set InSolveSet;
            public NSEquation[] Us;
            public NSEquation[] Vs;
            public NSEquation[] Ws;
            public SIMPLEPressure[] Ps;
            public SIMPLEVelocity[] Vels;

            protected int ZoneToSetRefP = -1;//<0则不需要设置
            public void CheckCreate(Set set, SIMPLEFields[] equations, EquationConnectBase[] connects)
            {
                this.InSolveSet = set;
                this.Equations = equations;
                this.Connects = connects;
                CheckCreator.ArrayStrict(ref Us, this.Equations.Length);
                CheckCreator.ArrayStrict(ref Vs, this.Equations.Length);
                CheckCreator.ArrayStrict(ref Ws, this.Equations.Length);
                CheckCreator.ArrayStrict(ref Ps, this.Equations.Length);
                CheckCreator.ArrayStrict(ref Vels, this.Equations.Length);
                for (int i1 = 0; i1 < this.Equations.Length; i1++)
                {
                    Us[i1] = this.Equations[i1].Velocity.U;
                    Vs[i1] = this.Equations[i1].Velocity.V;
                    Ws[i1] = this.Equations[i1].Velocity.W;
                    Ps[i1] = this.Equations[i1].Pressure;
                    Vels[i1] = this.Equations[i1].Velocity;
                }
                checkRefP();
            }
            protected void checkRefP()
            {
                bool isExistPressureBound()
                {
                    foreach (var v in Equations)
                    {
                        foreach (var b in v.SimpleBounds)
                        {
                            if (b is SIMPLEFields.PressureInOut)
                            {
                                return true;
                            }
                        }
                    }
                    return false;
                }
                if (isExistPressureBound())
                {
                    ZoneToSetRefP = -1;
                }
                else
                {
                    Vector3d centerCal = new Vector3d();
                    double totalVolumn = 0;
                    foreach (var p in this.Ps)
                    {
                        centerCal += p.Zone.Inner.Volumn * p.Zone.InnCells[p.ReferenceCell].Center;
                        totalVolumn += p.Zone.Inner.Volumn;
                    }
                    centerCal *= (1 / totalVolumn);
                    double minlenthSquare = double.MaxValue;
                    for (int i1 = 0; i1 < Ps.Count(); i1++)
                    {
                        ref var p = ref this.Ps[i1];
                        double ls = (p.Zone.InnCells[p.ReferenceCell].Center - centerCal).LengthSquared();
                        if (ls < minlenthSquare)
                        {
                            minlenthSquare = ls;
                            ZoneToSetRefP = i1;
                        }
                    }
                }
            }
            protected void setReferencePoint()
            {
                if (ZoneToSetRefP >= 0)
                {
                    Ps[ZoneToSetRefP].SetReferencePoint();
                }
            }
            protected void momentumEquation()
            {
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    foreach (var v in Vels)
                    {
                        v.FormAbOnFaces(threadInfo);
                    }
                });
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    foreach (var v in Vels)
                    {
                        v.FormAbOnCells(threadInfo);
                    }
                });
                switch (LxSim.Sim.Physics.TimeDiscreteType)
                {
                    case TimeDiscreteTypes.Steady: break;
                    case TimeDiscreteTypes.Implicit_1st:
                    case TimeDiscreteTypes.Implicit_2nd:
                        MultiThreads.RunParallelTask((threadInfo) =>
                        {
                            foreach (var v in Vels)
                            {
                                v.U.TransData.FormAb_ImplicitTrans(threadInfo);
                                v.V.TransData.FormAb_ImplicitTrans(threadInfo);
                                v.W.TransData.FormAb_ImplicitTrans(threadInfo);
                            }
                        });
                        break;
                    case TimeDiscreteTypes.Explicit_1st:
                    case TimeDiscreteTypes.Explicit_2nd:
                        MultiThreads.RunParallelTask((threadInfo) =>
                        {
                            foreach (var v in Vels)
                            {
                                v.U.TransData.FormAb_ExplicitTrans(threadInfo);
                                v.V.TransData.FormAb_ExplicitTrans(threadInfo);
                                v.W.TransData.FormAb_ExplicitTrans(threadInfo);
                            }
                        });
                        break;
                }
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    foreach (var v in Ps) { v.SaveDc(threadInfo); }
                });
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    foreach (var z in Vels)
                    {
                        z.U.ImplicitUnderRelax(threadInfo);
                        z.V.ImplicitUnderRelax(threadInfo);
                        z.W.ImplicitUnderRelax(threadInfo);
                    }
                });
                Solve(Us, Connects);
                Solve(Vs, Connects);
                Solve(Ws, Connects);
                Limit(Us, Connects);
                Limit(Vs, Connects);
                Limit(Ws, Connects);
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    foreach (var v in Equations) { v.Velocity.UpdateBoundsValue(threadInfo); }
                    foreach (var v in Equations) { v.MassFlux.CalFaceMassFlux(threadInfo); }
                });
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    foreach (var v in Equations) { v.MassFlux.CalMassImbalance(threadInfo); }
                });
            }
            protected void pressureEquation(bool isSimpleNotPiso)
            {
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    foreach (var pc in Ps) { pc.SaveXAndGradLastStep(threadInfo); }
                });
                for (int i1 = 0; i1 <= InSolveSet.NonOrthCorLoops; i1++)//即使不进行非正交修正，第一次（i1=0）也要运行
                {
                    if (i1 == 0)
                    {
                        MultiThreads.RunParallelTask((threadInfo) =>
                        {
                            foreach (var pc in Ps) { pc.FormABOnFaces(threadInfo); }
                        });
                        MultiThreads.RunParallelTask((threadInfo) =>
                        {
                            foreach (var pc in Ps) { pc.FormABOnCells(threadInfo); }
                        });
                    }
                    else
                    {
                        MultiThreads.RunParallelTask((threadInfo) =>
                        {
                            foreach (var pc in Ps) { pc.CalGradient_LSQ(threadInfo); }
                        });//如果是非正交修正，需要计算梯度；如果不进行非正交修正的话，后面的计算梯度和进行修正是同时进行的
                        MultiThreads.RunParallelTask((threadInfo) =>
                        {
                            foreach (var pc in Ps) { pc.UpdateABOnFaces(threadInfo); }
                        });
                        MultiThreads.RunParallelTask((threadInfo) =>
                        {
                            foreach (var pc in Ps) { pc.UpdateABOnCells(threadInfo); }
                        });
                    }
                    setReferencePoint();
                    Solve(Ps, Connects);
                    MultiThreads.RunParallelTask((threadInfo) =>
                    {
                        foreach (var pc in Ps) { pc.UpdateBoundsValue(threadInfo); }
                    });
                    foreach (var m in Ps) { m.LimitScalarValue(); }
                }
                if (isSimpleNotPiso)
                {
                    MultiThreads.RunParallelTask((threadInfo) =>
                    {
                        foreach (var pc in Ps) { pc.CalGradAndCorrect_SIMPLE_1(threadInfo); }
                    });
                    MultiThreads.RunParallelTask((threadInfo) =>
                    {
                        foreach (var pc in Ps) { pc.CalGradAndCorrect_SIMPLE_2(threadInfo); }
                    });
                }
                else
                {
                    MultiThreads.RunParallelTask((threadInfo) =>
                    {
                        foreach (var pc in Ps) { pc.CalGradAndCorrect_PISO(threadInfo); }
                    });
                }
                //更新速度边界、梯度等
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    foreach (var v in Equations) { v.Velocity.UpdateBoundsValue(threadInfo); }
                    foreach (var v in Equations) { v.MassFlux.CalFaceMassFlux(threadInfo); }
                });
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    foreach (var v in Equations) { v.MassFlux.CalMassImbalance(threadInfo); }
                });
            }
            public override void Initialize()
            {
                foreach (var e in Equations)
                {
                    e.Velocity.Initialize();
                    e.Pressure.Initialize();
                    e.MassFlux.Initialize();//massflux需要等pressure和velocity有了之后再计算
                }
            }
            public override void Solve_Steady()
            {
                for (int i2 = 0; i2 <= InSolveSet.AdditionalOuterLoops; i2++)
                {
                    momentumEquation();
                    pressureEquation(true);
                    for (int i1 = 0; i1 < InSolveSet.PISOLoops; i1++)
                    {
                        pressureEquation(false);
                    }
                }
                MultiThreads.RunParallelTask((threadInfo) =>
                {
                    foreach (var v in Equations)
                    {
                        //压力更新后也会影响到velocity，所以前面不更新了，最后更新一次
                        v.Velocity.U.CalGradient_LSQ(threadInfo);
                        v.Velocity.V.CalGradient_LSQ(threadInfo);
                        v.Velocity.W.CalGradient_LSQ(threadInfo);
                        //用于计算残差
                        v.MassFlux.CalResidual(threadInfo);
                        v.Velocity.U.CalResidual(threadInfo);
                        v.Velocity.V.CalResidual(threadInfo);
                        v.Velocity.W.CalResidual(threadInfo);
                    }
                });
            }
            public override void Solve_ImplicitTrans()
            {
                Solve_Steady();//momentumEquation中主动判断是否有时间
            }
            public override void SaveOldValue()
            {
                foreach (var v in Vels)
                {
                    v.U.TransData.SaveOldValue();
                    v.V.TransData.SaveOldValue();
                    v.W.TransData.SaveOldValue();
                }
            }
        }
    }
}
